function error_value = process_error_2d(mesh, pde, U, error_type)
gauss_weights_tri = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
psi = reference_basis_points_2d(mesh.basis_type);
dNdx = reference_basis_points_2d(mesh.basis_type, "dx");
dNdy = reference_basis_points_2d(mesh.basis_type, "dy");
ref_trans_tri = reference_basis_points_2d("P1");
switch error_type
    case "L2"
        sum_n = 0;
        for n = 1:mesh.Ne
            E = mesh.node(mesh.elem(n,:),:);
            V = E(1:3,:);
            uh_E = zeros(1,9);
            for j = 1:mesh.Nlb
                uh_E = uh_E + U(mesh.elem(n,j)).*psi{j};
            end
            ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
            f = (pde.u(ref_points(:,1),ref_points(:,2))' - uh_E).^2;
            r = (gauss_weights_tri * f') .* mesh.absJ;
            sum_n = sum_n + r;
        end
        error_value = sqrt(sum_n);
    case "H1"
        sum_nx = 0;
        sum_ny = 0;
        for n = 1:mesh.Ne
            E = mesh.node(mesh.elem(n,:),:);
            V = E(1:3,:);
            Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
            J_inv = inv(Jacobi);
            uh_E = zeros(2,9);
            phi = cellfun(@(x,y) J_inv' * [x; y], dNdx, dNdy, "UniformOutput", false);
            for j = 1:mesh.Nlb
                uh_E = uh_E + U(mesh.elem(n,j)) .* phi{j};
            end
            uh_Ex = uh_E(1,:);
            uh_Ey = uh_E(2,:);
            ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
            f1 = (pde.dudx(ref_points(:,1),ref_points(:,2))' - uh_Ex).^2;
            f2 = (pde.dudy(ref_points(:,1),ref_points(:,2))' - uh_Ey).^2;
            rx = (gauss_weights_tri * f1') .* mesh.absJ;
            ry = (gauss_weights_tri * f2') .* mesh.absJ;
            sum_nx = sum_nx + rx;
            sum_ny = sum_ny + ry;
        end
        error_value = sqrt(sum_nx + sum_ny);
    otherwise
        error("Invalid error type.");
end
end